####################################################################################
####################################################################################
########################BREEZE OR HURRICANE REPLICATION#############################
####################################################################################
####################################################################################

####################################################################################
########################STEP 1: PRELIMINARIES########################
####################################################################################

########1.1 DIRECTORY########
setwd("***DIRECTORY OF REPLICATION FOLDER***")

########1.2 LIBRARIES########
library(plyr)
library(dplyr)
library(zoo)
library(broom)
library(lubridate)
library(grid)
library(gridExtra)
library(foreign)
library(stargazer)
library(lattice)
library(lme4)
library(gsubfn)
library(sjPlot)
library(ggplot2)
library(interplot)
library("Hmisc")
library(panelAR)
library(dotwhisker)
library(broom)
library(dplyr)

########1.3 FUNCTIONS########
###creating shared legend for multiple plots###
grid_arrange_shared_legend <- function(..., ncol = length(list(...)), nrow = 1, position = c("bottom", "right")) {
  plots <- list(...)
  position <- match.arg(position)
  g <- ggplotGrob(plots[[1]] + 
                    theme(legend.position = position))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  lwidth <- sum(legend$width)
  gl <- lapply(plots, function(x) x +
                 theme(legend.position = "none"))
  gl <- c(gl, ncol = ncol, nrow = nrow)
  
  combined <- switch(position,
                     "bottom" = arrangeGrob(do.call(arrangeGrob, gl), 
                                            legend,ncol = 1,
                                            heights = unit.c(unit(1, "npc") - lheight, lheight)),
                     "right" = arrangeGrob(do.call(arrangeGrob, gl),
                                           legend, ncol = 2,
                                           widths = unit.c(unit(1, "npc") - lwidth, lwidth)))
  
  grid.newpage()
  grid.draw(combined)
  
  # return gtable invisibly
  invisible(combined)
}#function to arrange plots
###country cluster test formula
get_confint<-function(model, vcovCL){
  t<-qt(.975, model$df.residual)
  ct<-coeftest(model, vcovCL)
  est<-cbind(ct[,1], ct[,1]-t*ct[,2], ct[,1]+t*ct[,2])
  colnames(est)<-c("Estimate","LowerCI","UpperCI")
  return(est)
}
cluster.se<-function(model, cluster)
{
  require(multiwayvcov)
  require(lmtest)
  vcovCL<-cluster.vcov(model, cluster)
  coef<-coeftest(model, vcovCL)
  return(coef)
}


####################################################################################
########################STEP 2: DATA########################
####################################################################################

########2.1 READ IN DATA########
bh_rep_data <- read.csv("breeze_hurricane_replication_data.csv")#replication data for models
bh_rep_data$gwid <- as.factor(paste(bh_rep_data$gwid))#country code has to be factor for FE- models
populism_party_year <- read.csv("breeze_hurricane_populism_year.csv")#all populist parties in sample (yearly)
populism_party_periods <- read.csv("breeze_hurricane_populism_periods.csv")#all populist parties in sample (unique periods)
indicator <- read.csv("./db_indicator_description.csv",sep=",")

########2.2 Additional dataframes########
bh_rep_data_la <- subset(bh_rep_data, as.numeric(paste(gwid)) > 20 & as.numeric(paste(gwid)) < 200)#Latin America
bh_rep_data_eu <- subset(bh_rep_data, as.numeric(paste(gwid)) >= 200)#Europe
bh_rep_data_alt <- subset(bh_rep_data, perseat_lr_or_pop_altna_l1 / perseat_included_l1 < 0.2)#only country years where 80% of included party seat shares are coded as well with the more limited alternative measure

########2.3 DEFINE VARIABLES########
controls <- "regime_durability_l1 + gdppc_l1 + accession_perspective + post_crisis + year"#basic controls
controls2 <- "+ presidentialism_l1 + pr_elec_l1 + netmigration_perpop_l1 + oilrents_l1 + gini_l1"#expanded controls
principles <- c("FREEDOM","CONTROL","EQUALITY")
functions <- c("INDLIB", "RULEOFLAW", "PUBLIC", "COMPET", "MUTUCONS", "GOVCAP", "TRANSPAR", "PARTICIP", "REPRES")
components <- c("IL_PHIN", "IL_SELFU", "RL_EQL", "RL_QUAL", "PS_FRAS", "PS_FROP", "CO_COMP", "CO_OPEN", "MC_CHECKS", "MC_VERT", "GC_GORE", "GC_CEIM", "TR_NOSEC", "TR_PTPP", "PAR_EQPA", "PAR_EFPA", "REP_SR", "REP_DR")
subcomponents <- c("IL_PHIN1", "IL_PHIN2", "IL_PHIN3", "IL_SELFU1", "IL_SELFU2", "IL_SELFU3", "RL_EQL1", "RL_EQL2", "RL_EQL3", "RL_QUAL1", "RL_QUAL2", "RL_QUAL3", "PS_FRAS1", "PS_FRAS2", "PS_FRAS3", "PS_FROP1", "PS_FROP2", "PS_FROP3", "CO_COMP1", "CO_COMP2", "CO_COMP3", "CO_OPEN1", "CO_OPEN2", "CO_OPEN3", "MC_CHECKS1", "MC_CHECKS2", "MC_CHECKS3", "MC_VERT1", "MC_VERT2", "GC_GORE1", "GC_GORE2", "GC_GORE3", "GC_CEIM1", "GC_CEIM2", "GC_CEIM3", "GC_CEIM4", "TR_NOSEC1", "TR_NOSEC2", "TR_PTPP1", "TR_PTPP2", "TR_PTPP3", "PAR_EQPA1", "PAR_EQPA2", "PAR_EQPA3", "PAR_EFPA1", "PAR_EFPA2", "PAR_EFPA3", "REP_SR1", "REP_SR2", "REP_SR3", "REP_DR1", "REP_DR2", "REP_DR3")
indicators <- c("Consttort", "Convtort", "Politterr", "Torture", "Homicide", "Riot", "Constrel", "Constfreemov", "Freerelig", "Freemove", "Propright", "Secprop", "Constfair", "Pubtrial", "Judindepcor", "Judindepinf", "Impcourts", "Integrlegal", "Profjudg", "Proftenure", "Confjust", "Fairjust", "Confpolice", "Fairpolice", "Constfras", "Constass", "Union", "Memproorg", "Memhuman", "Memenviron", "Constspeech", "Constpress", "Newsimp", "Newspaper", "Balpress", "Neutrnp", "Meandistrict", "Gerryman", "Largpavo", "Votediff", "Herfindex", "Seatdiff", "Adminhurd", "Eff_thresh", "Smallpavo", "Enep", "Ceilings", "Funding", "Balexleg", "Balpowexle", "Seatsgov_2", "Judrev", "Powjudi", "Federalism", "Bicameralism", "Subexp", "Subrev", "Leglen", "Gov_term", "Confgov", "Devbehav", "Govstab", "Cabchange", "Antigovact", "Violantigov", "Mip", "Rip", "Publser", "Govdec", "Bureau", "CenBank_Ind", "Discinco", "Discexp", "Corrup", "CPI", "RestricFoi", "EffFoi", "Legmedia", "Polmedia", "Transp", "Suffrage", "Regprovap", "Repturnined", "Repturngeag", "Repaltined", "Repaltgeag", "Facilitat", "Registr", "Meanpart", "Eff_DD", "Petitions", "Demons", "Seatperinh", "No_districts", "Dirdem", "DD_Quora", "Gallagindex", "Issuecongr", "Polrightwom", "Constraints", "Partreg", "Womrep", "womgov", "Minrep", "Minpower")




####################################################################################
########################STEP 3: MAIN SPECIFICATIONS I + ROBUSTNESS CHECKS########################
####################################################################################

##########A: MAIN ARTICLE##########

########3.1 MODEL 1 [MAIN ARTICLE]########
models1 <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1 <- lm(", i, " ~ perseat_pop_l2_l1 + perseat_pop_c2_l1 + perseat_pop_r2_l1 + perseat_pop_gov_l1 + perseat_l2_l1 + perseat_c2_l1 + perseat_gov_l1 + ", controls, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1')) {models1[[",n,"]] <- ",i,"1} else {models1[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1)[n] <- i
  n <- n + 1
}

##########B: APPENDIX II: NESTED RESULTS##########

########3.2 MODEL 1a########
models1a <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1a <- lm(", i, " ~ perseat_pop_l1 + ", controls, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1a')) {models1a[[",n,"]] <- ",i,"1a} else {models1a[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1a)[n] <- i
  n <- n + 1
}

########3.3 MODEL 1b########
models1b <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1b <- lm(", i, " ~  perseat_pop_opp_l1 + perseat_pop_gov_l1 + perseat_gov_l1 + ", controls, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1b')) {models1b[[",n,"]] <- ",i,"1b} else {models1b[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1b)[n] <- i
  n <- n + 1
}

########3.4 MODEL 1c########
models1c <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1c <- lm(", i, " ~ perseat_pop_l2_l1 + perseat_pop_c2_l1 + perseat_pop_r2_l1 + perseat_l2_l1 + perseat_c2_l1 + ", controls, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1c')) {models1c[[",n,"]] <- ",i,"1c} else {models1c[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1c)[n] <- i
  n <- n + 1
}

##########C: APPENDIX IV: FURTHER ROBUSTNESS CHECKS##########

########3.5 MODEL 1 [Eur.]########
models1_eu <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1_eu <- lm(", i, " ~ perseat_pop_l2_l1 + perseat_pop_c2_l1 + perseat_pop_r2_l1 + perseat_pop_gov_l1 + perseat_l2_l1 + perseat_c2_l1 + perseat_gov_l1 + ", controls, " + gwid, data = bh_rep_data_eu)",sep="")
  step2 <- paste("if (exists('",i,"1_eu')) {models1_eu[[",n,"]] <- ",i,"1_eu} else {models1_eu[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1_eu)[n] <- i
  n <- n + 1
}

########3.6 MODEL 1a [Eur.]########
models1a_eu <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1a_eu <- lm(", i, " ~ perseat_pop_l1 + ", controls, " + gwid, data = bh_rep_data_eu)",sep="")
  step2 <- paste("if (exists('",i,"1a_eu')) {models1a_eu[[",n,"]] <- ",i,"1a_eu} else {models1a_eu[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1a_eu)[n] <- i
  n <- n + 1
}

########3.7 MODEL 1b [Eur.]########
models1b_eu <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1b_eu <- lm(", i, " ~  perseat_pop_opp_l1 + perseat_pop_gov_l1 + perseat_gov_l1 + ", controls, " + gwid, data = bh_rep_data_eu)",sep="")
  step2 <- paste("if (exists('",i,"1b_eu')) {models1b_eu[[",n,"]] <- ",i,"1b_eu} else {models1b_eu[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1b_eu)[n] <- i
  n <- n + 1
}

########3.8 MODEL 1c [Eur.]########
models1c_eu <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1c_eu <- lm(", i, " ~ perseat_pop_l2_l1 + perseat_pop_c2_l1 + perseat_pop_r2_l1 + perseat_l2_l1 + perseat_c2_l1 + ", controls, " + gwid, data = bh_rep_data_eu)",sep="")
  step2 <- paste("if (exists('",i,"1c_eu')) {models1c_eu[[",n,"]] <- ",i,"1c_eu} else {models1c_eu[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1c_eu)[n] <- i
  n <- n + 1
}

########3.9 MODEL 1a [LA.]########
models1a_la <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1a_la <- lm(", i, " ~ perseat_pop_l1 + ", controls, " + gwid, data = bh_rep_data_la)",sep="")
  step2 <- paste("if (exists('",i,"1a_la')) {models1a_la[[",n,"]] <- ",i,"1a_la} else {models1a_la[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1a_la)[n] <- i
  n <- n + 1
}

########3.10 MODEL 1 [alt. populism measure]########
models1_alt <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1_alt <- lm(", i, " ~ perseat_pop_alt_l2_l1 + perseat_pop_alt_c2_l1 + perseat_pop_alt_r2_l1 + perseat_pop_alt_gov_l1 + perseat_l2_l1 + perseat_c2_l1 + perseat_gov_l1 + ", controls, " + gwid, data = bh_rep_data_alt)",sep="")
  step2 <- paste("if (exists('",i,"1_alt')) {models1_alt[[",n,"]] <- ",i,"1_alt} else {models1_alt[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1_alt)[n] <- i
  n <- n + 1
}

########3.11 MODEL 1 [threshold 2]########
models1_t2 <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1_t2 <- lm(", i, " ~ perseat_pop_l3_l1 + perseat_pop_c3_l1 + perseat_pop_r3_l1 + perseat_pop_gov_l1 + perseat_l3_l1 + perseat_c3_l1 + perseat_gov_l1 + ", controls, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1_t2')) {models1_t2[[",n,"]] <- ",i,"1_t2} else {models4[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1_t2)[n] <- i
  n <- n + 1
}

########3.12 MODEL 1 [threshold 3]########
models1_t3 <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1_t3 <- lm(", i, " ~ perseat_pop_l4_l1 + perseat_pop_c4_l1 + perseat_pop_r4_l1 + perseat_pop_gov_l1 + perseat_l4_l1 + perseat_c4_l1 + perseat_gov_l1 + ", controls, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1_t3')) {models1_t3[[",n,"]] <- ",i,"1_t3} else {models4[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1_t3)[n] <- i
  n <- n + 1
}

########3.13 MODEL 1 [more controls]########
models1_controls <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1_controls <- lm(", i, " ~ perseat_pop_l2_l1 + perseat_pop_c2_l1 + perseat_pop_r2_l1 + perseat_pop_gov_l1 + perseat_l2_l1 + perseat_c2_l1 + perseat_gov_l1 + ", controls, controls2, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1_controls')) {models1_controls[[",n,"]] <- ",i,"1_controls} else {models1_controls[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1_controls)[n] <- i
  n <- n + 1
}

########3.14 MODEL 1 [separate ideology terms gov./opp.]########
models1_inter <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1_inter <- lm(", i, " ~ perseat_pop_opp_l2_l1 + perseat_pop_gov_l2_l1 + perseat_pop_opp_c2_l1 + perseat_pop_gov_c2_l1 + perseat_pop_opp_r2_l1 + perseat_pop_gov_r2_l1 + perseat_l2_l1 + perseat_c2_l1 + perseat_gov_l1 + ", controls, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1_inter')) {models1_inter[[",n,"]] <- ",i,"1_inter} else {models4[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1_inter)[n] <- i
  n <- n + 1
}

########3.15 MODEL 1 [Prais-Winsten]########
models1_prais <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"1_prais <- panelAR(", i, " ~ perseat_pop_l2_l1 + perseat_pop_c2_l1 + perseat_pop_r2_l1 + perseat_pop_gov_l1 + perseat_l2_l1 + perseat_c2_l1 + perseat_gov_l1 + ", controls, " + gwid, data = bh_rep_data, panelVar = 'gwid', timeVar = 'year', autoCorr='ar1', rho.na.rm=TRUE, panel.weight='t-1', bound.rho=TRUE,panelCorrMethod='pcse')",sep="")
  step2 <- paste("if (exists('",i,"1_prais')) {models1_prais[[",n,"]] <- ",i,"1_prais} else {models1_prais[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1_prais)[n] <- i
  n <- n + 1
}


####################################################################################
########################STEP 4: MAIN SPECIFICATIONS II + ROBUSTNESS CHECKS########################
####################################################################################

##########A: MAIN ARTICLE##########

########4.1 MODEL 2 [MAIN ARTICLE / APPENDIX III]########
models2 <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"2 <- lm(", i, " ~ perseat_pop_largest_opp_l1 + perseat_pop_largest_gov_l1 + lr_pop_largest_l1 + I(lr_pop_largest_l1^2)  + perseat_largest_opp_l1 + perseat_largest_gov_l1 +  lr_largest_l1 + I(lr_largest_l1^2) + ", controls, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1')) {models2[[",n,"]] <- ",i,"2} else {models2[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models2)[n] <- i
  n <- n + 1
}

##########B: APPENDIX III: NESTED RESULTS##########

########4.2 MODEL 2a########
models2a <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"2a <- lm(", i, " ~ perseat_pop_largest_l1 + perseat_largest_l1 + ", controls, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1')) {models2a[[",n,"]] <- ",i,"2a} else {models2a[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models2a)[n] <- i
  n <- n + 1
}

########5.3 MODEL 2b########
models2b <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"2b <- lm(", i, " ~ perseat_pop_largest_opp_l1 + perseat_pop_largest_gov_l1 + perseat_largest_opp_l1 + perseat_largest_gov_l1 + ", controls, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1')) {models2b[[",n,"]] <- ",i,"2b} else {models2b[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models2b)[n] <- i
  n <- n + 1
}

########5.4 MODEL 2c########
models2c <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"2c <- lm(", i, " ~ perseat_pop_largest_l1 + lr_pop_largest_l1 + I(lr_pop_largest_l1^2) + perseat_largest_l1 + lr_largest_l1 + I(lr_largest_l1^2) + ", controls, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1')) {models2c[[",n,"]] <- ",i,"2c} else {models2c[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models2c)[n] <- i
  n <- n + 1
}

##########C: APPENDIX IV: FURTHER ROBUSTNESS CHECKS##########

########5.5 MODEL 2 [Cubic Terms]########
models2_cub <- list()
n <- 1
for (i in (c(functions))) {
  step1 <- paste(i,"2_cub <- lm(", i, " ~ perseat_pop_largest_opp_l1 + perseat_pop_largest_gov_l1 + lr_pop_largest_l1 + I(lr_pop_largest_l1^2) + I(lr_pop_largest_l1^3) + perseat_largest_opp_l1 + perseat_largest_gov_l1 +  lr_largest_l1 + I(lr_largest_l1^2) + I(lr_largest_l1^3) + ", controls, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1')) {models2_cub[[",n,"]] <- ",i,"2_cub} else {models2_cub[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models2_cub)[n] <- i
  n <- n + 1
}


####################################################################################
########################STEP 5: EXPORT ALL REGRESSION RESULTS (MAIN + APPENDIX)########################
####################################################################################

##########A: MAIN ARTICLE##########

########5.1 MODEL 1 [MAIN ARTICLE]########
stargazer(COMPET1, PARTICIP1, REPRES1, GOVCAP1, INDLIB1, RULEOFLAW1, MUTUCONS1, TRANSPAR1, PUBLIC1, covariate.labels = c("Populist seat share (left)","Populist seat share (centre)","Populist seat share (right)","Populist seat share (gov.)","Seat share (left)","Seat share (centre)","Seat share (gov.)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table 3: Model 1 - Populist party strength (left-centre-right), government inclusion, and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")

########5.2 MODEL 2 [MAIN ARTICLE / APPENDIX III]########
stargazer(COMPET2, PARTICIP2, REPRES2, GOVCAP2, INDLIB2, RULEOFLAW2, MUTUCONS2, TRANSPAR2, PUBLIC2, covariate.labels = c("Populist seat share (largest, opp.)","Populist seat share (largest, gov.)","Populist left-right (largest)","Populist left-right (largest) sq.","Seat share (largest, opp.)","Seat share (largest, gov.)","Left-right (largest)","Left-right (largest) sq.", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A7: Model 2 - Largest populist party strength (gov. vs. opp., left-right position) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")

##########B: APPENDIX##########

########5.3 APPENDIX II: NESTED MODEL RESULTS########
stargazer(COMPET1a, PARTICIP1a, REPRES1a, GOVCAP1a, INDLIB1a, RULEOFLAW1a, MUTUCONS1a, TRANSPAR1a, PUBLIC1a, covariate.labels = c("Populist seat share", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A4: Model 1a - Populist party strength and the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET1b, PARTICIP1b, REPRES1b, GOVCAP1b, INDLIB1b, RULEOFLAW1b, MUTUCONS1b, TRANSPAR1b, PUBLIC1b, covariate.labels = c("Populist seat share (opp.)","Populist seat share (gov.)","Seat share (gov.)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A5: Model 1b - Populist party strength (gov. vs. opp.) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET1c, PARTICIP1c, REPRES1c, GOVCAP1c, INDLIB1c, RULEOFLAW1c, MUTUCONS1c, TRANSPAR1c, PUBLIC1c, covariate.labels = c("Populist seat share (left)","Populist seat share (centre)","Populist seat share (right)","Seat share (left)","Seat share (centre)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A6: Model 1c - Populist party strength (left-centre-right) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")

########5.4 APPENDIX III: ALTERNATIVE SPECIFICATION (NESTED RESULTS)########
stargazer(COMPET2a, PARTICIP2a, REPRES2a, GOVCAP2a, INDLIB2a, RULEOFLAW2a, MUTUCONS2a, TRANSPAR2a, PUBLIC2a, covariate.labels = c("Populist seat share (largest)","Seat share (largest)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A8: Model 2a  - Strongest populist party strength and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET2b, PARTICIP2b, REPRES2b, GOVCAP2b, INDLIB2b, RULEOFLAW2b, MUTUCONS2b, TRANSPAR2b, PUBLIC2b, covariate.labels = c("Populist seat share (largest, opp.)","Populist seat share (largest, gov.)","Seat share (largest, opp.)","Seat share (largest, gov.)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A9: Model 2b - Largest populist party strength (gov. vs. opp.) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET2c, PARTICIP2c, REPRES2c, GOVCAP2c, INDLIB2c, RULEOFLAW2c, MUTUCONS2c, TRANSPAR2c, PUBLIC2c, covariate.labels = c("Populist seat share (largest)","Populist left-right (largest)","Populist left-right (largest) sq.","Seat share (largest)","Left-right (largest)","Left-right (largest) sq.", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A10: Model 2c - Largest populist party strength (left-right position) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")

########5.5 APPENDIX IV: FURTHER ROBUSTNESS CHECKS########
stargazer(COMPET1_eu, PARTICIP1_eu, REPRES1_eu, GOVCAP1_eu, INDLIB1_eu, RULEOFLAW1_eu, MUTUCONS1_eu, TRANSPAR1_eu, PUBLIC1_eu, covariate.labels = c("Populist seat share (left)","Populist seat share (centre)","Populist seat share (right)","Populist seat share (gov.)","Seat share (left)","Seat share (centre)","Seat share (gov.)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A11: Model 1 [Eur.] - Populist party strength (gov. vs. opp., left-centre-right) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET1a_eu, PARTICIP1a_eu, REPRES1a_eu, GOVCAP1a_eu, INDLIB1a_eu, RULEOFLAW1a_eu, MUTUCONS1a_eu, TRANSPAR1a_eu, PUBLIC1a_eu, covariate.labels = c("Populist seat share", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A12: Model 1a [Eur.] - Populist party strength and the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET1b_eu, PARTICIP1b_eu, REPRES1b_eu, GOVCAP1b_eu, INDLIB1b_eu, RULEOFLAW1b_eu, MUTUCONS1b_eu, TRANSPAR1b_eu, PUBLIC1b_eu, covariate.labels = c("Populist seat share (opp.)","Populist seat share (gov.)","Seat share (gov.)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A13: Model 1b [Eur.] - Populist party strength (gov. vs. opp.) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET1c_eu, PARTICIP1c_eu, REPRES1c_eu, GOVCAP1c_eu, INDLIB1c_eu, RULEOFLAW1c_eu, MUTUCONS1c_eu, TRANSPAR1c_eu, PUBLIC1c_eu, covariate.labels = c("Populist seat share (left)","Populist seat share (centre)","Populist seat share (right)","Seat share (left)","Seat share (centre)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A14: Model 1c [Eur.] - Populist party strength (left-centre-right) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET1a_la, PARTICIP1a_la, REPRES1a_la, GOVCAP1a_la, INDLIB1a_la, RULEOFLAW1a_la, MUTUCONS1a_la, TRANSPAR1a_la, PUBLIC1a_la, covariate.labels = c("Populist seat share", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A15: Model 1a [LA.] - Populist party strength and the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET1_alt, PARTICIP1_alt, REPRES1_alt, GOVCAP1_alt, INDLIB1_alt, RULEOFLAW1_alt, MUTUCONS1_alt, TRANSPAR1_alt, PUBLIC1_alt, covariate.labels = c("Populist seat share (left)","Populist seat share (centre)","Populist seat share (right)","Populist seat share (gov.)","Seat share (left)","Seat share (centre)","Seat share (gov.)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A16: Model 1 [alt. populism measure] - Populist party strength (gov. vs. opp., left-centre-right) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET1_t2, PARTICIP1_t2, REPRES1_t2, GOVCAP1_t2, INDLIB1_t2, RULEOFLAW1_t2, MUTUCONS1_t2, TRANSPAR1_t2, PUBLIC1_t2, covariate.labels = c("Populist seat share (left)","Populist seat share (centre)","Populist seat share (right)","Populist seat share (gov.)","Seat share (left)","Seat share (centre)","Seat share (gov.)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A17: Model 1 [threshold 2] - Populist party strength (gov. vs. opp., left-centre-right) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET1_t3, PARTICIP1_t3, REPRES1_t3, GOVCAP1_t3, INDLIB1_t3, RULEOFLAW1_t3, MUTUCONS1_t3, TRANSPAR1_t3, PUBLIC1_t3, covariate.labels = c("Populist seat share (left)","Populist seat share (centre)","Populist seat share (right)","Populist seat share (gov.)","Seat share (left)","Seat share (centre)","Seat share (gov.)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A18: Model 1 [threshold 3] - Populist party strength (gov. vs. opp., left-centre-right) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET1_controls, PARTICIP1_controls, REPRES1_controls, GOVCAP1_controls, INDLIB1_controls, RULEOFLAW1_controls, MUTUCONS1_controls, TRANSPAR1_controls, PUBLIC1_controls, covariate.labels = c("Populist seat share (left)","Populist seat share (centre)","Populist seat share (right)","Populist seat share (gov.)","Seat share (left)","Seat share (centre)","Seat share (gov.)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year","Presidentialism","PR","Net migration","Oil exports","GINI", "Constant"), title = "Table A19: Model 1 [more controls] - Populist party strength (gov. vs. opp., left-centre-right) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET1_inter, PARTICIP1_inter, REPRES1_inter, GOVCAP1_inter, INDLIB1_inter, RULEOFLAW1_inter, MUTUCONS1_inter, TRANSPAR1_inter, PUBLIC1_inter, covariate.labels = c("Populist seat share (left, opp.)","Populist seat share (left, gov.)","Populist seat share (centre, opp.)","Populist seat share (centre, gov.)","Populist seat share (right, opp.)","Populist seat share (right, gov.)","Seat share (left)","Seat share (centre)","Seat share (gov.)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A20: Model 1 [separate ideology terms gov./opp.] - Populist party strength (gov. vs. opp., left-centre-right) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
stargazer(COMPET2_cub, PARTICIP2_cub, REPRES2_cub, GOVCAP2_cub, INDLIB2_cub, RULEOFLAW2_cub, MUTUCONS2_cub, TRANSPAR2_cub, PUBLIC2_cub, covariate.labels = c("Populist seat share (largest, opp.)","Populist seat share (largest, gov.)","Populist left-right (largest)","Populist left-right (largest) sq.","Populist left-right (largest) cub.","Seat share (largest, opp.)","Seat share (largest, gov.)","Left-right (largest)","Left-right (largest) sq.","Left-right (largest)","Left-right (largest) cub.", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A21: Model 2 [Cubic Term] - Largest populist party strength (gov. vs. opp., left-right position) and components of the quality of democracy.", colnames = F, dep.var.labels =c("Competition","Participation","Representation", "Government capability", "Individual liberties", "Rule of law","Constraints","Transparency","Public sphere"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")
#Table A22 (COMPET1_prais, PARTICIP1_prais, REPRES1_prais, GOVCAP1_prais, INDLIB1_prais, RULEOFLAW1_prais, MUTUCONS1_prais, TRANSPAR1_prais, PUBLIC1_prais) manually created


####################################################################################
########################STEP 6: TABLES + FIGURES (MAIN PART)########################
####################################################################################

########6.1 Table II [populism periods]########
tableII <- subset(bh_rep_data, !is.na(perseat_lr_or_popna_l1))
tableII <- unique(tableII[c("country","gwid")])#all countries in our analysis frame
tableII$gwid <- as.integer(paste(tableII$gwid))
tableII <- left_join(tableII,populism_party_year[c("gwid","pname","mrcode","year","perseat","govparty","populism","lr")], by=c("gwid"))
tableII_lag <- tableII
tableII_lag[c("country","pname","perseat","govparty","lr")] <- NULL
colnames(tableII_lag)[4] <- paste(colnames(tableII_lag)[4], "l1", sep = "_")
tableII_lag$year <- tableII_lag$year + 1
tableII <- left_join(tableII, tableII_lag, by =c("gwid","mrcode","year"))
tableII$lr <- round(tableII$lr,2)
tableII$perseat <- round(tableII$perseat,2)
tableII$period <- ifelse(is.na(tableII$populism_l1) | tableII$populism != tableII$populism_l1, 1, 0)
tableII <- tableII %>% group_by(gwid,mrcode) %>% arrange(gwid,mrcode,year) %>% mutate(period = cumsum(period))#period number increases everytime the populism dichotomous value changes
tableII <- tableII %>% group_by(gwid,mrcode, period) %>% arrange(gwid,mrcode,year) %>% mutate(from = min(year))#year in which the period starts
tableII <- tableII %>% group_by(gwid,mrcode, period) %>% arrange(gwid,mrcode,year) %>% mutate(to = max(year))#year in which the period ends
tableII <- tableII %>% group_by(gwid,mrcode, period) %>% arrange(gwid,mrcode,year) %>% mutate(perseat_max = max(perseat, na.rm=T))#maximum seat shares in period
tableII <- tableII %>% group_by(gwid,mrcode, period) %>% arrange(gwid,mrcode,year) %>% mutate(gov_max = max(govparty, na.rm=T))#govparty in period
tableII <- unique(tableII[c("country","pname","from", "to","perseat_max","gov_max","populism","lr")])
tableII <- tableII[order(-as.numeric(tableII$perseat_max)),] 
tableII <- tableII[order(-as.numeric(tableII$populism)),] 
tableII <- unique(tableII)
write.csv(tableII,file="./figures/main/tableII.csv")

########6.2 Figure I [scatterplot]########
bh_rep_data_p_plot <- subset(populism_party_periods,populism == 1)
figureI <- ggplot(bh_rep_data_p_plot, aes(x=lr, y=perseat,colour=factor(govparty))) +  geom_point(alpha = 0.5) + geom_segment(aes(x=-0.25,xend=-0.25,y= 0,yend=1), colour = 'black', linetype="dashed") + geom_segment(aes(x=0.25,xend=0.25,y= 0,yend=1), colour = 'black', linetype="dashed")+ scale_color_manual(values = c('#fd8d3c',"#b10026","black"),labels=c("Populist in opp.","Populist in gov."),name="Party type") + xlab('Left-right stance') + ylab('Seat share') + theme_bw()  +theme(text=element_text(family="Times"))#+ ggtitle("Distribution of cases by populist seat shares, ideology, and government inclusion")
ggsave(file="./figures/main/figureI.jpg", figureI, width = 17, height = 6.5, units ="cm", dpi = 800)

########6.3 Figure II: underlying data [map]########
figureII <- populism_party_periods[c("country","gwid","lr","populism","perseat")]
figureII$rightpop <- ifelse(figureII$populism == 1 & figureII$lr >= 0.25,1,0)
figureII$leftpop <- ifelse(figureII$populism == 1 & figureII$lr <= -0.25,1,0)
figureII$centrepop <- ifelse(figureII$populism == 1 & figureII$lr > -0.25 & figureII$lr < 0.25,1,0)
figureII$lr_available <- ifelse(!is.na(figureII$lr), 1, 0)
figureII <- figureII %>% group_by(gwid) %>% mutate(perseat_maxpop = max(populism * perseat * figureII$lr_available,na.rm=T))
figureII <- figureII %>% group_by(gwid) %>% mutate(populism_ever = max(populism,na.rm=T))
figureII <- figureII %>% group_by(gwid) %>% mutate(populism_everavailable = max(populism * lr_available,na.rm=T))
figureII <- subset(figureII, perseat * populism * lr_available == perseat_maxpop | (populism_ever == 0) | (populism_ever == 1 & populism_everavailable == 0))
figureII$perseat <- ifelse(figureII$populism_ever == 0, 0, figureII$perseat)
figureII$rightpop <- ifelse(figureII$populism_ever == 0, 0, figureII$rightpop)
figureII$leftpop <- ifelse(figureII$populism_ever == 0, 0, figureII$leftpop)
figureII$centrepop <- ifelse(figureII$populism_ever == 0, 0, figureII$centrepop)
figureII <- unique(figureII[c("country","gwid","rightpop","leftpop","centrepop","perseat","populism_ever","populism_everavailable")])
figureII_r <- subset(figureII, rightpop == 1)
figureII_l <- subset(figureII, leftpop == 1)
figureII_c <- subset(figureII, centrepop == 1)
figureII_none <- subset(figureII, populism_ever == 0)
write.csv(figureII_r,"./figures/main/figureII/figureII_r.csv")
write.csv(figureII_l,"./figures/main/figureII/figureII_l.csv")
write.csv(figureII_c,"./figures/main/figureII/figureII_c.csv")
write.csv(figureII_none,"./figures/main/figureII/figureII_none.csv")

########6.4 Figure III [coefficient plot]########
for (i in (c(functions))) {
  step1 <- paste(i,"1_df <- data.frame(subset(tidy(",i,"1, conf.int = T, conf.level = 0.9), term == 'perseat_pop_l2_l1' | term == 'perseat_pop_c2_l1' | term == 'perseat_pop_r2_l1' | term == 'perseat_pop_gov_l1'))",sep="")
  step2 <- paste(i,"1_df$model2 <- ",i,"1_df$term",sep="")
  eval(parse(text=c(step1, step2)))
}
COMPET1_df$term <- "Competition"
PARTICIP1_df$term <- "Participation"
REPRES1_df$term <- "Representation"
GOVCAP1_df$term <- "Government capability"
INDLIB1_df$term <- "Ind. liberties"
RULEOFLAW1_df$term <- "Rule of law"
MUTUCONS1_df$term <- "Constraints"
TRANSPAR1_df$term <- "Transparency"
PUBLIC1_df$term <- "Public sphere"
figureIII <- rbind(COMPET1_df,PARTICIP1_df,REPRES1_df,GOVCAP1_df,INDLIB1_df,RULEOFLAW1_df,MUTUCONS1_df,TRANSPAR1_df,PUBLIC1_df)
figureIII$model <- ifelse(figureIII$model2 == "perseat_pop_l2_l1", "Left pop.", NA)
figureIII$model <- ifelse(figureIII$model2 == "perseat_pop_c2_l1", "Center pop.", figureIII$model)
figureIII$model <- ifelse(figureIII$model2 == "perseat_pop_r2_l1", "Right pop.", figureIII$model)
figureIII$model <- ifelse(figureIII$model2 == "perseat_pop_gov_l1", "Gov. pop.", figureIII$model)
figureIII$model <- as.factor(figureIII$model)
figureIII$model = factor(figureIII$model,levels(figureIII$model)[c(2,4,1,3)])
figureIII$model = factor(figureIII$model,levels(figureIII$model)[c(4,3,2,1)])
figureIII <- figureIII%>% arrange(desc(model))
figureIII_plot <- dwplot(figureIII,
                       vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) +
  theme_bw() + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("") +
  scale_color_manual(name = "Coefficient for:", values = c("black","#0571b0","#fec44f","#d7191c"),
                     breaks=c("Left pop.","Center pop.","Right pop.","Gov. pop."))+ theme(legend.position = "bottom", text=element_text(family="Times"))
ggsave(file="./figures/main/figureIII.jpg", figureIII_plot, width = 20.1, height = 15, units ="cm", dpi = 800)

########6.5 Figure IV [predictions plot]########
testdata_nopop <- data.frame(perseat_pop_largest_opp_l1 = 0, perseat_pop_largest_gov_l1 = 0, lr_pop_largest_l1 = 0, perseat_largest_opp_l1 = 0.05359, perseat_largest_gov_l1 = 0.3651, lr_largest_l1 = 0.12759, regime_durability_l1 = 38.76, gdppc_l1 = 24640, oilrents_l1 = 0.83704, accession_perspective = 0, year = 2016, post_crisis = 1, gwid = '200')
testdata_nongov <- data.frame(perseat_pop_largest_opp_l1 = rep(0.5, each = 800), perseat_pop_largest_gov_l1 = 0, lr_pop_largest_l1 =  rep(seq(-1,1,by=0.005012531), 2), perseat_largest_opp_l1 = 0.05359, perseat_largest_gov_l1 = 0.3651, lr_largest_l1 = 0.12759, regime_durability_l1 = 38.76, gdppc_l1 = 24640, oilrents_l1 = 0.83704, accession_perspective = 0, year = 2016, post_crisis = 1, gwid = '200')
testdata_nongov$popingov <- 0
testdata_gov <- data.frame(perseat_pop_largest_opp_l1 = 0, perseat_pop_largest_gov_l1  = rep(0.5, each = 800), lr_pop_largest_l1 =  rep(seq(-1,1,by=0.005012531), 2), perseat_largest_opp_l1 = 0.05359, perseat_largest_gov_l1 = 0.3651, lr_largest_l1 = 0.12759, regime_durability_l1 = 38.76, gdppc_l1 = 24640, oilrents_l1 = 0.83704, accession_perspective = 0, year = 2016, post_crisis = 1, gwid = '200')
testdata_gov$popingov <- 1
testdata_all <- rbind(testdata_nongov, testdata_gov)
testdata_all$perseat_pop_largest_either <- pmax(testdata_all$perseat_pop_largest_opp_l1, testdata_all$perseat_pop_largest_gov_l1)
for (i in (c(functions))) {
  step1 <- paste("test_nopop_",i," <- predict(",i,"2, newdata = testdata_nopop, interval = 'confidence',level =0.9)",sep="")#predict situation if no populists
  step2 <- paste("testdata_nopop_full_",i," <- cbind(testdata_nopop, test_nopop_",i,")",sep="")#add to df
  eval(parse(text=c(step1,step2)))
  step3 <- paste("test_nongov_",i," <- predict(",i,"2, newdata = testdata_all, interval = 'confidence',level =0.9)",sep="")#predict situation if populists, situations as above
  step4 <- paste("testdata_nongov_full_",i," <- cbind(testdata_all, test_nongov_",i,")",sep="")
  step5 <- paste(i, "_figureIV <- ggplot(data = testdata_nongov_full_",i,") + xlim(-1,1) + geom_line(aes(x=lr_pop_largest_l1, y=fit, linetype = factor(popingov), colour = factor(perseat_pop_largest_either)),size=0.5) + xlab('Populist ideology (left-right)') + ylab('Predicted value') + ggtitle ('",i,"')  + scale_color_manual(values=c('#fc4e2a','#e31a1c','#b10026'), name ='Populist seat shares')   + scale_linetype_manual(values=c('longdash','solid','dotted'), labels = c('not in gov.','in gov.','none'), name ='Populists')  +
                 geom_ribbon(data = subset(testdata_nongov_full_",i,", popingov == 0), aes(x=lr_pop_largest_l1, ymin = lwr, ymax = upr), fill = 'grey', alpha = 0.2) +
                 geom_ribbon(data = subset(testdata_nongov_full_",i,", popingov == 1), aes(x=lr_pop_largest_l1, ymin = lwr, ymax = upr), fill = 'grey', alpha = 0.2) +
                 
                 scale_size_manual( values = c(2,4,2,1) ) + theme_bw() +theme(text=element_text(family='Times')) +
                 geom_segment(aes(x=-1,xend=1,y= testdata_nopop_full_",i,"$fit,yend=testdata_nopop_full_",i,"$fit,linetype = 'dotted'))",sep="")
  eval(parse(text=c(step3, step4, step5)))
}
figureIV <- grid_arrange_shared_legend(COMPET_figureIV + ggtitle("Competition"), PARTICIP_figureIV + ggtitle("Participation"), REPRES_figureIV + ggtitle("Representation"), GOVCAP_figureIV + ggtitle("Government capability"), INDLIB_figureIV + ggtitle("Individual liberties"), RULEOFLAW_figureIV + ggtitle("Rule of law"), MUTUCONS_figureIV + ggtitle("Mutual constraints"), TRANSPAR_figureIV + ggtitle("Transparency"), PUBLIC_figureIV + ggtitle("Public sphere"), ncol=3,nrow=3)
ggsave(file="./figures/main/figureIV.jpg", figureIV, width = 20.1, height = 13.8, units ="cm", dpi = 800)


####################################################################################
########################STEP 7: TABLES + FIGURES (APPENDIX)########################
####################################################################################

##########A: APPENDIX I##########

########7.1 Table A1 [Correlations of populism measures]########
populism_correlations <- populism_party_periods[c("mrcode","gwid","populism1","populism2","populism3","populism4","populism5","populism6","populism7","populism8")]
for (i in (c("populism1","populism2","populism3","populism4","populism5","populism6","populism7","populism8"))) {
  step1 <- paste("populism_correlations <- populism_correlations %>% group_by(gwid,mrcode) %>% arrange(gwid,mrcode) %>% mutate(",i,"=max(",i,",na.rm=T))",sep="")
  eval(parse(text=c(step1)))#collects list of (numeric) variables
}
is.na(populism_correlations) <- sapply(populism_correlations, is.infinite)
populism_correlations <- unique(populism_correlations[c("mrcode","gwid","populism1","populism2","populism3","populism4","populism5","populism6","populism7","populism8")])
populism_correlations[c("mrcode","gwid")] <- NULL
tablea1 <- rcorr(as.matrix(populism_correlations))

########7.1 Table A3 [Correlations of ideology measures]########
lr_correlations <- unique(populism_party_periods[c("mrcode","gwid","lr_cm","lr_hi","lr_bl","lr_chess","lr_bg","lr_manuscript_transformed")])
lr_correlations[c("mrcode","gwid")] <- NULL
tablea3 <- rcorr(as.matrix(lr_correlations))

##########B: APPENDIX IV##########

########7.3 Figure A1 [Predictions of model 2_cubic]########
testdata_nopop <- data.frame(perseat_pop_largest_opp_l1 = 0, perseat_pop_largest_gov_l1 = 0, lr_pop_largest_l1 = 0, perseat_largest_opp_l1 = 0.05359, perseat_largest_gov_l1 = 0.3651, lr_largest_l1 = 0.12759, regime_durability_l1 = 38.76, gdppc_l1 = 24640, oilrents_l1 = 0.83704, accession_perspective = 0, year = 2016, post_crisis = 1, gwid = '200')
testdata_nongov <- data.frame(perseat_pop_largest_opp_l1 = rep(0.5, each = 800), perseat_pop_largest_gov_l1 = 0, lr_pop_largest_l1 =  rep(seq(-1,1,by=0.005012531), 2), perseat_largest_opp_l1 = 0.05359, perseat_largest_gov_l1 = 0.3651, lr_largest_l1 = 0.12759, regime_durability_l1 = 38.76, gdppc_l1 = 24640, oilrents_l1 = 0.83704, accession_perspective = 0, year = 2016, post_crisis = 1, gwid = '200')
testdata_nongov$popingov <- 0
testdata_gov <- data.frame(perseat_pop_largest_opp_l1 = 0, perseat_pop_largest_gov_l1  = rep(0.5, each = 800), lr_pop_largest_l1 =  rep(seq(-1,1,by=0.005012531), 2), perseat_largest_opp_l1 = 0.05359, perseat_largest_gov_l1 = 0.3651, lr_largest_l1 = 0.12759, regime_durability_l1 = 38.76, gdppc_l1 = 24640, oilrents_l1 = 0.83704, accession_perspective = 0, year = 2016, post_crisis = 1, gwid = '200')
testdata_gov$popingov <- 1
testdata_all <- rbind(testdata_nongov, testdata_gov)
testdata_all$perseat_pop_largest_either <- pmax(testdata_all$perseat_pop_largest_opp_l1, testdata_all$perseat_pop_largest_gov_l1)
for (i in (c(functions))) {
  step1 <- paste("test_nopop_",i," <- predict(",i,"2_cub, newdata = testdata_nopop, interval = 'confidence',level =0.9)",sep="")#predict situation if no populists
  step2 <- paste("testdata_nopop_full_",i," <- cbind(testdata_nopop, test_nopop_",i,")",sep="")#add to df
  eval(parse(text=c(step1,step2)))
  step3 <- paste("test_nongov_",i," <- predict(",i,"2_cub, newdata = testdata_all, interval = 'confidence',level =0.9)",sep="")#predict situation if populists, situations as above
  step4 <- paste("testdata_nongov_full_",i," <- cbind(testdata_all, test_nongov_",i,")",sep="")
  step5 <- paste(i, "_figurea1 <- ggplot(data = testdata_nongov_full_",i,") + xlim(-1,1) + geom_line(aes(x=lr_pop_largest_l1, y=fit, linetype = factor(popingov), colour = factor(perseat_pop_largest_either)),size=0.5) + xlab('Populist ideology (left-right)') + ylab('Predicted value') + ggtitle ('",i,"')  + scale_color_manual(values=c('#fc4e2a','#e31a1c','#b10026'), name ='Populist seat shares')   + scale_linetype_manual(values=c('longdash','solid','dotted'), labels = c('not in gov.','in gov.','none'), name ='Populists')  +
                 geom_ribbon(data = subset(testdata_nongov_full_",i,", popingov == 0), aes(x=lr_pop_largest_l1, ymin = lwr, ymax = upr), fill = 'grey', alpha = 0.2) +
                 geom_ribbon(data = subset(testdata_nongov_full_",i,", popingov == 1), aes(x=lr_pop_largest_l1, ymin = lwr, ymax = upr), fill = 'grey', alpha = 0.2) +
                 
                 scale_size_manual( values = c(8,4,2,1) ) + theme_bw() +theme(text=element_text(family='Times')) +
                 geom_segment(aes(x=-1,xend=1,y= testdata_nopop_full_",i,"$fit,yend=testdata_nopop_full_",i,"$fit,linetype = 'dotted'))",sep="")
  eval(parse(text=c(step3, step4, step5)))
}
figurea1 <- grid_arrange_shared_legend(COMPET_figurea1 + ggtitle("Competition"), PARTICIP_figurea1 + ggtitle("Participation"), REPRES_figurea1 + ggtitle("Representation"), GOVCAP_figurea1 + ggtitle("Government capability"), INDLIB_figurea1 + ggtitle("Individual liberties"), RULEOFLAW_figurea1 + ggtitle("Rule of law"), MUTUCONS_figurea1 + ggtitle("Mutual constraints"), TRANSPAR_figurea1 + ggtitle("Transparency"), PUBLIC_figurea1 + ggtitle("Public sphere"), ncol=3,nrow=3)
ggsave(file="./figures/appendix/figurea1.jpg", figurea1, width = 20.1, height = 13.8, units ="cm", dpi = 800)

##########B: APPENDIX V##########

########7.3 Table A23 [placebo tests]########
perseat_pop_p <- lm(perseat_pop ~  COMPET_l1 + PARTICIP_l1 + REPRES_l1 + PUBLIC_l1 + INDLIB_l1 + RULEOFLAW_l1 + MUTUCONS_l1 + TRANSPAR_l1 + GOVCAP_l1 + regime_durability_l1 + gdppc_l1 + accession_perspective + post_crisis + year + gwid, data = bh_rep_data)
perseat_pop_l2_p <- lm(perseat_pop_l2 ~  COMPET_l1 + PARTICIP_l1 + REPRES_l1 + PUBLIC_l1 + INDLIB_l1 + RULEOFLAW_l1 + MUTUCONS_l1 + TRANSPAR_l1 + GOVCAP_l1 + perseat_l2_l1 + perseat_c2_l1 + regime_durability_l1 + gdppc_l1 + accession_perspective + post_crisis + year + gwid, data = bh_rep_data)
perseat_pop_c2_p <- lm(perseat_pop_c2 ~  COMPET_l1 + PARTICIP_l1 + REPRES_l1 + PUBLIC_l1 + INDLIB_l1 + RULEOFLAW_l1 + MUTUCONS_l1 + TRANSPAR_l1 + GOVCAP_l1 +  perseat_l2_l1 + perseat_c2_l1 + regime_durability_l1 + gdppc_l1 + accession_perspective + post_crisis + year + gwid, data = bh_rep_data)
perseat_pop_r2_p <- lm(perseat_pop_r2 ~  COMPET_l1 + PARTICIP_l1 + REPRES_l1 + PUBLIC_l1 + INDLIB_l1 + RULEOFLAW_l1 + MUTUCONS_l1 + TRANSPAR_l1 + GOVCAP_l1 + perseat_l2_l1 + perseat_c2_l1 + regime_durability_l1 + gdppc_l1 + accession_perspective + post_crisis + year + gwid, data = bh_rep_data)
stargazer(perseat_pop_p, perseat_pop_l2_p, perseat_pop_c2_p, perseat_pop_r2_p, covariate.labels = c("Competition","Participation","Representation","Public sphere","Individual liberties","Rule of law","Mutucal constraints","Transparency","Government capability","Seat shares (left)","Seat shares (centre)", "Democratic consolidation", "GDP p.c.", "Credible EU accession perspective", "Post-crisis", "Year", "Constant"), title = "Table A22: Placebo tests", colnames = F, dep.var.labels =c("Populist seat share", "Populist seat share (left)", "Populist seat share (centre)", "Populist seat share (right)"), omit=c("gwid"), add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes","Yes")), type = "text")

##########C: APPENDIX VI##########

########7.4 Table A24 [indicator-level results]########
models1 <- list()
n <- 1
for (i in (c(indicators))) {
  step1 <- paste(i,"1 <- lm(", i, " ~ perseat_pop_l2_l1 + perseat_pop_c2_l1 + perseat_pop_r2_l1 + perseat_pop_gov_l1 + perseat_l2_l1 + perseat_c2_l1 + perseat_gov_l1 + ", controls, " + gwid, data = bh_rep_data)",sep="")
  step2 <- paste("if (exists('",i,"1')) {models1[[",n,"]] <- ",i,"1} else {models1[[",n,"]] <- NA}",sep="")
  eval(parse(text=c(step1, step2)))
  names(models1)[n] <- i
  n <- n + 1
}
models1 <- ldply(models1, tidy, .id = "model")
tablea24 <- subset(models1, term == "perseat_pop_l2_l1" | term == "perseat_pop_c2_l1" | term == "perseat_pop_r2_l1" | term == "perseat_pop_gov_l1")
tablea24 <- left_join(tablea24,indicator[c("model","short_description","function.","subfunction","component")])
tablea24 <- tablea24[order(tablea24$function.,tablea24$subfunction,tablea24$component),]
tablea24$star <- ifelse(tablea24$p.value < 0.1, paste("*"), "")
tablea24$star <- ifelse(tablea24$p.value < 0.05, paste("**"), tablea24$star)
tablea24$star <- ifelse(tablea24$p.value < 0.01, paste("***"), tablea24$star)
tablea24$summary <- paste(round(tablea24$estimate, digits = 2), " (",round(tablea24$std.error, digits = 2),")",tablea24$star,sep="")
tablea24$term_long <- ifelse(tablea24$term == "perseat_pop_l2_l1", "Populist seat share (left)", NA)
tablea24$term_long <- ifelse(tablea24$term == "perseat_pop_c2_l1", "Populist seat share (centre)", tablea24$term_long)
tablea24$term_long <- ifelse(tablea24$term == "perseat_pop_r2_l1", "Populist seat share (right)", tablea24$term_long)
tablea24$term_long <- ifelse(tablea24$term == "perseat_pop_gov_l1", "Populist seat share (gov.)", tablea24$term_long)
tablea24 <- tablea24[c("model","short_description","term_long","summary","function.")]
tablea24 <- reshape(tablea24, idvar = c("model","short_description","function."), timevar = "term_long", direction = "wide")
write.csv(tablea24,file="./figures/appendix/tablea24.csv")